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Abstract 

We calculate the inclusive decay rate of the spin-triplet bottomonium states Xbj into charm 
hadrons, including the leading-order color-singlet and color-octet bb annihilation mechanisms. We 
also calculate the momentum distribution of the charm quark from the decay of Xbj- The infrared 
divergences from the color-singlet process bb — ► ccg are factored into the probability density at 
the origin for a bb pair in a color-octet state. That probability density can be determined phe- 
nomenologically from the fraction of decays of Xbj that include charm hadrons. It can then be 
used to predict the partial widths into light hadrons for all four states in the P-wave bottomonium 
multiplet. 

PACS numbers: 12.38.-t, 12.39.St, 13.20. Gd, 14.40. Gx 
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I. INTRODUCTION 



The asymptotic freedom of QCD suggests that the total widths of heavy quarkonium 
states should be calculable using perturbation theory. The earliest calculations of the widths 
of P-wave quarkonium states using perturbative QCD were plagued with infrared divergences 
[1-3]. The calculations were based on a factorization assumption that the width could be 
expressed as the product of |P'(0)| 2 , where R'(0) is the derivative of the radial wave function 
at the origin, and a perturbatively calculable coefficient. However the coefficients were found 
to be infrared divergent at leading order in a s for the spin-1 states and at next-to-leading 
order in a s for the spin-0 and spin-2 states. The infrared divergences were often expressed 
in terms of a logarithmic dependence on the binding energy of the quarkonium, a quantity 
that is not calculable using perturbation theory. However the correct interpretation of the 
infrared divergences is that they reveal the failure of the factorization assumption. 

This problem was overcome in 1992 when Bodwin, Braaten, and Lepage showed that the 
infrared divergences could be absorbed into the probability for the heavy-quark-antiquark 
(QQ) P a i r to be at the same point in a color-octet state [4]. They used a nonrelativistic 
effective field theory for the QQ sector of QCD called NRQCD to derive a general factor- 
ization formula for inclusive quarkonium decay rates [5]. A P-wave multiplet consists of 
four heavy quarkonium states: Xqo, Xq\i XQ2, and hq with J PC quantum numbers ++ , 
1 ++ , 2 ++ , and l + ~, respectively. At leading order in the velocity v of the heavy quark or 
antiquark in the quarkonium rest frame, there are only two independent nonperturbative 
factors in the annihilation decay rates of all four states in the P-wave multiplet: (Ci), which 
is proportional to |P'(0)| 2 , and (O s ), which is proportional to the probability for the Q and 
Q to be at the same point in a color-octet state. These nonperturbative factors can be 
expressed as matrix elements of local four-quark operators in NRQCD. The short-distance 
coefficients of the NRQCD matrix elements can be calculated as power series in the QCD 
coupling constant a s . 

The widths of all four states in a P-wave multiplet can be calculated by using the NRQCD 
factorization formula, once the two nonperturbative factors (Oi) and (C 8 ) have been deter- 
mined. These matrix elements can be calculated by using lattice simulations of NRQCD. An 
alternative is to estimate the color-singlet matrix element (0\) by using potential models 
and to determine the color-octet matrix element (Os) phenomenologically. The phenomeno- 



2 



logical determination of (0%) requires the measurement of an observable that is sensitive to 
this matrix element. In the case of bottomonium, one such observable is the inclusive rate 
for charm production in decays of the spin-triplet P-wave states Xbj- This rate is sensitive 
to (Og) because the production of charm quarks from bb annihilation in the color-singlet 
channel is suppressed by a factor of a s , relative to production in the color-octet channel. 

There has been little previous work on open charm production in bottomonium decays. 
In 1978, Fritzsch and Streng calculated the decay rate of T into charm at leading order 
in a s [6]. In 1979, Barbieri, Caffo, and Remiddi calculated the decay rates of the P-wave 
bottomonium states into charm at leading order in a s under the assumption that the rates 
could be expressed as products of |P'(0)| 2 and a perturbatively calculable coefficient [7]. In 
the case of xm an d Xb2, the coefficients contained infrared divergences that were expressed 
in terms of logarithms of the binding energy. However, as we have mentioned, the correct 
interpretation of the infrared divergences is that they are contained in the probability to find 
the QQ pair at a point in a color-octet state. By making use of the NRQCD factorization 
formalism, one can now carry out rigorous calculations of inclusive charm production from 
Xbj decays. 

In their 1979 paper, Barbieri, Caffo, and Remiddi calculated the invariant mass distribu- 
tion of the cc pair in Xbj decays. In order to make contact with experiment, one might be 
tempted to identify this distribution with the invariant mass distribution of pairs of charm 
hadrons. However that distribution cannot be measured easily because the probability of 
identifying both charm hadrons is very low. Furthermore, the effects of the hadronization of 
the charm quark into a charm hadron have a large effect on the distribution. These effects 
cannot be calculated perturbatively, and they would also be very difficult to measure. A 
more useful quantity to calculate is the momentum distribution of the charm quark in Xbj 
decays. This cannot be compared directly with the momentum distribution of the charm 
hadrons because of the effects of hadronization. However, the effects of hadronization can 
be determined experimentally by measuring the momentum distribution of charm hadrons 
in e + e~ annihilation. 

On the experimental side, the spin-triplet members of two multiplets of P-wave bottomo- 
nium states have been discovered: Xw(l-P) an d Xw(2P)- The only properties of these states 
that have been measured thus far are their masses and their radiative branching fractions 
into the S-wave bottomonium states T(nS). The total widths of the Xbj(nP) states have not 
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been measured. Recent runs of the CLEO experiment at the T(2S') and T(3S) resonances 
have provided new data on the Xbj(lP) and Xw(2P) states. The P-factory experiments 
BABAR and Belle can study the Xbj(nP) states by using data samples of T(2S) and T(3S) 
provided by initial-state radiation. The Belle experiment has also accumulated data by 
running directly on the T(3S) state. 

In this paper, we study inclusive charm production in P-wave bottomonium decays. In 
Sec. II, we present the NRQCD factorization formulas for the annihilation decays of P-wave 
bottomonium states, and we discuss the NRQCD matrix elements that appear as long- 
distance factors in the factorization formulas. In Sec. Ill, we calculate the charm-quark 
momentum distribution in decays of the spin-triplet P-wave states Xbj- We include the 
color-singlet process bb — > ccg, which has a short-distance coefficient of order o^, and the 
color-octet process bb — > cc, which has a short- distance coefficient of order a 2 . In Sec. IV, 
we calculate the inclusive rate into charm by integrating over the charm-quark momentum 
distribution. In Sec. V, we illustrate the momentum distribution for a charm meson D 
by convolving the charm-quark momentum distribution with a fragmentation function for 
c — > D that has been measured in e + e~ annihilation. Details of the calculations are presented 
in appendices. 

II. ANNIHILATION DECAYS OF P-WAVE BOTTOMONIUM 

The NRQCD factorization formula expresses the annihilation contribution to the hadronic 
width of a heavy quarkonium state as an infinite sum of products of short-distance coef- 
ficients, which can be calculated as power series in a s , and nonperturbative long-distance 
factors [5]. The long-distance factors can be expressed as expectation values of local four- 
quark operators O c ( 2S+1 Lj) that are defined in Ref. [5]. These NRQCD matrix elements 
scale as definite powers of the velocity v of the heavy quark in the quarkonium rest frame. 
For each of the P-wave states, there are only two matrix elements that contribute up to 
corrections of relative order v 2 : (Oi( 3 Pj)) XbJ and {0 8 ( 3 Si)) XbJ for Xbj and (Oi(}Pi)) hb and 
(Osi 1 S ))h b for hb- Heavy-quark spin symmetry can be used to reduce all these matrix el- 
ements at leading order in v to two independent matrix elements that we will denote by 



4 



{0,) Xb and (0 8 )$>: 

{Oi) Xb = (O^Pi))^ « {O^Pj)) Xb , v (la) 
(Os)^ = WSotfg « (0 8 ( 3 5i)>gJ. (lb) 

The superscript (A) on (0 8 )^ indicates the sensitivity of this matrix element to the NRQCD 
factorization scale. There is a total of 10 independent matrix elements that contribute 
through order v 2 [8]. 

The NRQCD factorization formulas for the annihilation widths of the Xbj at leading order 
in v can be expressed as 

T[ XbJ - *] = Aj(A) <^ + ^ (2) 

ml mi 

where X represents all possible states that consist of hadrons lighter than the B meson, and 
A is the NRQCD factorization scale. An analogous equation holds for the rate dT[xbj — > -X"] 
that is differential in the kinematic variables. The short-distance coefficients whose leading 
terms are order a 2 are 

a, - fg™2. < 3b > 

A 8 = ^rifira 2 , (3c) 

where N c = 3 is the number of colors, Cp = (N 2 — 1)/ (2A^ C ) = 4/3, ny = 4 is the number 
of light flavors of quarks, including charm, and the masses of the light quarks have been 
neglected. The coefficients A and A 2 were first calculated by Barbieri, Gatto, and Kogerler 
in 1976 [1]. The coefficient A 8 was first calculated for massless quarks in Ref. [5]. The 
short-distance coefficients whose leading terms are order are 



A 1 (A) = 



C F a 
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^" 288" ) CA+ {-Y7-9 l0g 2m-J n r 
where Ca = N c = 3, and, again, the masses of the quarks, including the charm quark, 
have been neglected. The coefficient A 1 was calculated in Refs. [9, 10]. The coefficients 
Aj(A) depend on A, beginning at order a^, in such a way as to cancel the dependence of 
the matrix element {Og)xV on A. The next-to-leading-order terms in the coefficients A , A 2 , 
and A 8 have been calculated by Petrelli, Cacciari, Greco, Maltoni, and Mangano [9] and by 



Huang and Chao [10]. The short-distance coefficients are insensitive to m c , the mass of the 
charm quark. The dependence of the leading terms in Aj and A$ on m c will be calculated 
in Sec. IV. The leading correction term in A 8 is proportional to a 2 s {m c /m b ) A . The leading 
correction terms in Aj are proportional to a^ s {m c /m} } ) 2 . 

In the NRQCD factorization formula in Eq. (2), the decay rates are summed over all light 
hadronic states. In most cases, there are no factorization formulas for less inclusive decay 
rates. An exception is the inclusive charm decay rate. The decay of Xbj into a final state 
that includes charm hadrons requires the annihilation of the bb pair into partons that include 
a cc pair. The mass of the charm quark is large enough that the contribution to the short- 
distance coefficients from bb annihilation into cc pairs may be calculable in perturbation 
theory. At leading order in v, the NRQCD factorization formula for the inclusive charm 
decay rate of Xbj involves the same matrix elements as the completely inclusive annihilation 
decay rate in Eq. (2): 



Y[ XbJ - c + X] = Af(A) + Af> (5) 



(Oi) Xb Ac) (0 8 )^ 
ml mi 



where c + X represents all possible states that include a charm hadron. The short-distance 
l^ c) and Af 



coefficients A^ and A^ are power series in a s whose coefficients are functions of the mass 



ratio m c /mt). We can deduce the limit as m c — > of the leading term in A^ from the value 
of A 8 in Eq. (3c): A^ — > (1/3)7^. Unlike the coefficients in the fully inclusive factorization 
formula in Eq. (2), the coefficients A^ and A^ in Eq. (5) are sensitive to the charm-quark 
mass. The leading terms in A^ and A^ will be calculated in Sec. IV. We will find that the 
leading term in A^j\ which is of order a^, depends logarithmically on m c /mb- 

The NRQCD matrix elements in Eqs. (1) can, in principle, be calculated by using lattice 
simulations of NRQCD. The feasibility of such calculations was first demonstrated by Bod- 
win, Sinclair, and Kim using quenched lattice NRQCD [11]. The best calculations available 
to date have been carried out using two dynamical light quarks [12]. After extrapolation to 
three light-quark flavors [12], the values for the IP multiplet are 

(Oi) Xb (iP) = 3.2±0.7GeV 5 , (6a) 
(O ) (A) 

X6(1P) - 0.0021 ± 0.0007 GeV~ 2 . (6b) 



We have estimated the errors for the three-flavor case by treating the systematic errors from 
the quenched and two-flavor calculations as 100% correlated, treating the statistical errors as 
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uncorrelated, and adding the resulting systematic and statistical errors for the three-flavor 
case in quadrature. The matrix element (OsJ^nw m Eq. (6b) was computed at A =4.3 
GeV. 

The color-singlet matrix elements can also be estimated by using potential models for 
heavy quarkonium: 

<^iW)^Kp(0)| 2 , (7) 

where N c = 3 is the number of colors and R n p(r) is the radial wave function for the nP 
multiplet. The values of \R' n p(0)\ 2 for four potential models have been tabulated in Ref. [13]. 
Using the value of |i?^ P (0)| 2 for the Buchmuller-Tye potential, we obtain 

(Oi) Xbi iP) » 2.03 GeV 5 , (8a) 
(0i)x*(2P) » 2.37 GeV 5 . (8b) 

The values of (Ci) Xi ,( n p) from the four potential models in Ref. [13] range from those in 
Eqs. (8) to those for the Cornell potential, which are about 50% larger. In the case of S- 
wave states, there has been recent progress in determining the color-singlet NRQCD matrix 
element from potential models [14]. The value of the radial wavefunction at the origin 
|-Ris(0)| 2 of the T(IS') that follows from these methods agrees most closely with that from 
the Buchmuller-Tye potential. 

With the choice of normalization of the operators in Ref. [5], the color-octet matrix 
element (0 8 ) Xb can be interpreted intuitively as the probability density at the origin for the 
bb pair to be in a color-octet state. One can obtain an order-of-magnitude estimate of a lower 
bound on the quantity (O s ) Xb by using the renormalization properties of the operators [5]. 
The operator O s depends on a renormalization scale A, and it mixes under renormalization 
with 0\. The solution to the renormalization group equation at leading order in a s is [5] 

where (3q = (lliV c — 2ny)/6 = 25/6 is the first coefficient in the beta function for QCD with 
71/ = 4 flavors of light quarks. The second term on the right side of Eq. (9) has a physical 
interpretation as a correction from gluon radiation that arises from gluon energies between 
A and m^. Since m^v is the typical momentum scale in a quarkonium state, we choose 
A = m b v. Unless there is a near cancellation between the two terms in Eq. (9) for A = m b v, 
the matrix element (Os)^^ should either be comparable to or larger than the second term 
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on the right side. This gives us an order-of-magnitude estimate of a lower bound on the 
matrix element: 

(0 8 ><"*> > — log ( i^ik (10) 

Since the one-loop bottom-quark pole mass is m^ olc ^ ps 4.6 GeV, we set mj = 4.6 GeV and 
rribV = 1.5 GeV. Then our estimated lower bound on the dimensionless ratio of the matrix 
elements 

Ps=ml(0 8 )^/(0 1 ) Xb (11) 

is ps > 0.068. In comparison, the lattice results in Eq. (6), taken with = 4.6 GeV, give 
Ps = 0.044 ± 0.015. Given the errors, this result is compatible with the estimated lower 
bound from Eq. (10). 

In NRQCD, there is no general relation between the matrix elements (Oi) Xb and (C^xT^ 
for different P-wave multiplets. However, if the scale m b v 2 is below the QCD scale A QCD , 
then the ratio p 8 in Eq. (11) is the same for all the P-wave multiplets [15]. 

III. CHARM QUARK PRODUCTION IN Xb DECAY 
A. Pert ur bat ive matching 

The coefficients in the NRQCD factorization formula for inclusive charm production in 
Eq. (5) are short- distance quantities that are insensitive to the long-distance behavior of 
the external bb states. This implies that the short-distance coefficients can be computed 
in perturbation theory. It also implies that, for purposes of computing the short- distance 
coefficients, we can replace the external bb hadronic states in the factorization formula with 
perturbative bb states. We compute the short-distance coefficients by matching the perturba- 
tive expressions for the bb annihilation rates in full QCD with the corresponding perturbative 
NRQCD factorization expressions for the annihilation rates. The perturbative analog of the 
NRQCD factorization formula in Eq. (5) for the annihilation rates of appropriate bb states 
is 

dT[bb ^c + X] = j2 dAf{K) (0l(3 y )fe5 + dA^ WS ^\ (12) 

J=0 m b m b 

We have written the factorization formula in differential form so that we can consider dis- 
tributions in kinematic variables associated with the charm quark. We can determine the 



four short-distance coefficients dA^ and dA^ by (i) calculating the annihilation rate in 
perturbative QCD for a bb pair in four appropriate independent bb states, (ii) calculating 
the NRQCD matrix elements for each of those four states using perturbative NRQCD, and 
then (iii) solving the linear set of equations for the coefficients. 

We wish to calculate the short-distance coefficients at leading order in a s , which is order 
a 2 for A^ and order oP s for A^f . At this order, we must take into account the renormalization 
of the NRQCD matrix element {O s ( 3 Si)) b t. We regularize the NRQCD matrix element 
by using dimensional regularization in d = 4 — 2e space-time dimensions, and we define 
the renormalized NRQCD matrix element by using the modified minimal subtraction (MS) 
prescription. The relation between the bare operator Og( 3 Si) and the renormalized operator 
8 ( 3 Si) (A) with NRQCD factorization scale A is [5, 9] 

8 (* Sl ) = O s CS^ + (4 " e ~ 7) \ 2 ^ 2 £ 0^ Pj ) + .... (13) 

e UV 3nN c m 2 j-jj 

The subscript UV indicates that the pole in e is associated with an ultraviolet divergence. 
We have shown explicitly only those terms that contribute through order a s and at leading 
order in v to the expectation values in a color-singlet P-wave bb state or in a color-octet 
S-wave bb state. 

The perturbative matrix elements of the NRQCD operators regularized with dimensional 
regularization are particularly simple if we also use dimensional regularization to regularize 
infrared divergences and we expand the matrix elements in powers of the relative momentum 
q of the b and b. In this case, all loop corrections to the regulated matrix element vanish 
because there is no scale for the dimensionally regularized integrals. In particular, the 
ultraviolet poles in e cancel the infrared poles in e. Thus, we have 



(OsCS,))^ = (0 8 ( 3 Si)>£ ee) > (14) 

where (0 8 ( 3 S'i))^ cg ' ) is the matrix element of the bare NRQCD operator with both infrared 
and ultraviolet divergences dimensionally regulated, and (C 8 ( 3 S'i))^ ree - ) is the tree-level ap- 
proximation to the matrix element of the bare NRQCD operator. If we take the expectation 
value of Eq. (13) in a bb state, dimensionally regulating both UV and IR divergences, and 
substitute (14), we find that 

= - (4 " R T> ' 3 2 ^ E< '' 3 ^)>^ + • • • • < 15 » 



The subscript IR indicates that the pole in e is now associated with an infrared divergence. 

To determine the four short-distance coefficients and in Eq. (5), we must calculate 
the annihilation rate for four appropriate bb states using perturbative QCD. A convenient 
choice for these states consists of a bb pair in a color-octet 3 Si state, which we denote by 
bbs( 3 Si), and a bb pair in each of the three color-singlet 3 Pj states, which we denote by 
bbii^Pj). For these states, the factorization formula in Eq. (12) reduces at leading order in 
a s and at leading order in v to 

dTibb^S,) ^c + X] = dAf ( ° a(3gl)) ^Wi) (16a) 

drfb^Pj) ^ c + x] = dAf(A) gig^jW^) + dA f (i 6b) 

m\ mi 

In Eq. (16a), the term involving the color-singlet operator does not contribute because it is 
of higher order in a s . In Eq. (16b), the color-octet matrix element (Os^Si)}^^^ can be 
simplified by using the fact that the tree-level term in Eq. (15) does not contribute. The 
factorization formula in Eq. (16b) can then be reduced to 

driftf J>,) - c + x] = f iAf W - (±L£Z>!^i iA f\ <^('ft)>*c-p,) , (17) 

V e m ott1\ c ) m b 

Eqs. (16a) and (17) can be solved to obtain the short- distance coefficients dA^ and dA^f 
in terms of the perturbative decay rates dT[bb 8 ( 3 Si) — > c + X] and dT[bbi ( 3 Pj) — > c + X] 
and the perturbative matrix elements (Os^Si)}^^^ and (Oi^Pj^^^Pj). At the order 
in a s of the present calculation, the perturbative matrix elements can be computed at tree 
level. In the next three subsections, we compute the required perturbative decay rates and 
perturbative matrix elements. As we will see, <ir[66 1 ( 3 Pj) — > c + X] contains an infrared 
divergence that is canceled by the explicit infrared divergence in the second term on the 
right side of Eq. (17). The short-distance coefficients are then infrared finite, as expected. 



B. Amplitudes for 66 annihilation into charm 

The momenta of the b and b that annihilate to produce charm can be expressed as 

V = \P + q, (18a) 
p = \P-q, (18b) 
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where P and q are the total and relative momenta of the bb pair. In the rest frame of the bb 
pair, the explicit momenta are P = (2E b , 0) and q = (0, q), where E b = y/m% + q 2 and m b 
is the mass of the bottom quark. An annihilation amplitude can be expressed in the form 

v(p)Au(p) = Tr[Au(p)v(p)], (19) 

where A is a matrix that acts on spinors with both Dirac and color indices. The amplitude 
in Eq. (19) can be projected into a particular spin and color channel by replacing u(p)v(p) 
with a projection matrix. The color projectors ~K\ and 7Tg onto a color-singlet state and onto 
a color-octet state with color index a are 

TTi = ^=1, (20a) 



< = V2T a , (20b) 

where 1 is the 3x3 unit matrix and T a is a generator of the fundamental representation 
of SU(3). The color projectors are normalized so that Tr[7r 1 7r{] = l and Tr[7Tg7rg ] = 5 ab . The 
projector onto a spin-triplet state with four-momentum P M , rest energy \fP 2 ~ = 2E b , and 
spin polarization vector es satisfying P ■ es = is es^IIg [16-18], where 

^= - ■ -(p , + m b )(P+2E b ) 1 »(P-m b ). (21) 

4V2E b (E b + m b ) 

The spin projector is normalized so that 

Tr[( e5 -n 3 )(e 5 -n 3 ) t ] =4p Po- (22) 

At leading order in v, the amplitude for the annihilation of a bb pair in a color-octet 
spin-triplet S-wave state with spin polarization vector es is es^Ag 1 , where 

A? = Tr[A(I%®*S)] • (23) 

The leading color-octet mechanism for producing charm in bb annihilation is via the process 
bb — > cc, whose rate is of order a 2 s . The matrix A for the process b{p)b{p) — > c(pi)c(p 2 ) is 

A[bb - cc] = ^( Pl )T b 7^(p 2 ) [T 6 7 ^] . (24) 

Using Eq. (23), we find that the coefficient of in the annihilation amplitude is 

AT = ^u(Pi)T a Yv(p 2 ), (25) 
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where we have omitted terms proportional to P M because P ■ es = 0. 

At leading order in the relative velocity v of the b or b in the quarkonium rest frame, 
the amplitude for the annihilation of a bb pair in a color-singlet spin-triplet P-wave state 
with spin polarization vector es and orbital-angular- momentum polarization vector is 
clvCs^" , where 

A F = jfc Tr l A ( n s®*i)]U- ( 26 ) 
The leading color-singlet mechanism for producing charm in bb annihilation is the process 
bb — > ccg, whose rate is of order a 3 . The matrix A for the process b(p)b(p) — > c(pi)c(p 2 )(7(p 3 ) 
is 

A[bb^ccg] = ^—^u( Pl )T a lx v(P2)e b ;( P 3) 



X 



[T a TV A(p - p 3 ) 7 * + T b T a ^K{-p + p 3 )7 A ] , (27) 



where A(k) is defined by 



Using Eq. (26), we find that 



A(*) = (28) 



AT = 2 ^J_ p3)2 HPi)T a ^v(P2)eT(Ps) 

x ^-Tr{ [ 7 A A(p - p 3 )Y + rM-P + Ps)l X ] ^}\ q=0 - (29) 

C. Color-octet short-distance coefficient 

(c) 

We proceed to calculate the differential coefficient dA 8 of the color-octet term in the 
NRQCD factorization formula. We use the perturbative factorization formula in Eq. (16a), 
which requires calculating the annihilation rate of a bb pair in a color-octet 3 Si state. The 
resulting expression for dA^ will also be needed in the determination of the coefficients 
dAj^ that makes use of the perturbative factorization formula in Eq. (17). In that equation, 
dA^ is multiplied by a pole in e. It is therefore necessary to calculate dA^ in d = 4 — 2e 
space-time dimensions. 

The differential annihilation rate of a color-octet 3 Si bb state into charm through the 
color-octet process bb — > cc can be expressed in the form 

dr[bb s ( 3 S 1 ) ^c + X]= (^L_l m ^ATAr^j d$ 2 , (30) 
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where A^ 1 is the amplitude in Eq. (25), d$ 2 is the differential 2-body phase space for cc, 
and is the projection tensor for spin 1: 

pH pv 

I" v = -<T + -JPT- (31) 

The factor of l/(d — 1) in Eq. (30) comes from averaging over the spin states of the bb pair. 
The explicit sum in Eq. (30) is over the color and spin states of the c and c. The evaluation 
of that sum gives 

I, a AT AT* = (47ra s A 2e ) 2 (iV c 2 - 1) (d - 2 + r) , (32) 

cc 

where g 2 s = 4na s A 2e and A is the scale associated with dimensional regularization. In our 
calculation, A becomes the NRQCD factorization scale. We have set E b — * m b in for 
consistency with the prescription for A^ in Eq. (23), which involved expanding to leading 
order in v. 

We wish to obtain an expression for the coefficient that is differential in the energy of 
the charm quark. We therefore integrate over the entire 2-body phase space, except for E 1: 
the energy of the charm quark in the bb rest frame. In the center-of-momentum frame, the 
differential 2-body phase space in d = 4 — 2e space-time dimensions reduces to 

d$2 = c 2 (e)^A(£i - E b )dE u (33) 

where = (Ef — m 2 ) 1 / 2 is the magnitude of the three- momentum of the charm quark and 
2E b is the energy of the bb pair. The dimensionless coefficient c 2 (e), which reduces to 1 as 
e — > 0, is defined by 



r 



-3^ 



It is useful to express the differential phase space in terms of an energy fraction x\ for the 
charm quark defined by 

Xl = E 1 /E b . (35) 

There is some ambiguity in the choice of E b . The choice E b = M XbJ /2 gives the correct 
kinematic limits on the energy of the charm quark. However, we choose E b = m b in order 
to maintain consistency with the nonrelativistic approximation that we used in computing 
A^ 1 in Eq. (23). The expression for the differential phase space then reduces to 



[(1 - r)miY 8n 
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where r is the square of the ratio of the charm- and bottom-quark masses: 

r = ml/ml. (37) 

Inserting the differential phase space in Eq. (36) into Eq. (30) and using Eq. (32), we find 
that the expression for the differential annihilation rate reduces to 

dT 66 8 ( 3 S'i) -> c + X = ; x ^— - — — ^ (1-xi 38 

[(1 — r)m£\ e ct — 1 

To complete the matching calculation of the coefficient dA^\ we need to evaluate the 
NRQCD matrix element on the right side of the factorization formula in Eq. (16a). The bb 
states have the standard relativistic normalizations. At leading order in the nonrelativistic 
expansion, the matrix element is therefore 

(08( 3 5i)) fe 5 8( 35 1 ) = 4(iV c 2 -lK. (39) 

Inserting Eqs. (38) and (39) into Eq. (16a), we find that the differential coefficient dA^ in 
d dimensions is 

, Ac) c 2 (e)A 4e (d-2 + r)Vl -r 2 ... N , ,. n . 
= [(l-rKl» X 2(d-l) < 4() ) 

Upon setting e = 0, we find that the differential coefficient with respect to the energy fraction 

of the charm quark reduces to 

jg = < i+ y^ (i-,,). 

D. Color-singlet short-distance coefficients 

We next calculate the differential coefficients dAy of the color-singlet terms in the 
NRQCD factorization formula. We use the perturbative factorization formula in Eq. (17), 
which requires calculating the annihilation rate of a bb pair in a color-singlet 3 Pj state for 
J = 0, 1, and 2. This annihilation rate is infrared divergent at leading order in a s . We use 
dimensional regularization in d = 4 — 2e space-time dimensions to regularize the infrared 
divergence. 

The differential annihilation rate of a color-singlet 3 Pj bb state into charm through the 
color-singlet process bb — > ccg can be expressed in the form 

dribb^Pj) ^c + x] = ^L^K^rf E atat^j dt> 3 , (42) 
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where A^ 1 is the amplitude in Eq. (29), rf$ 3 is the differential three-body phase space for ccg, 
and the K J a ^ are the projection tensors for total angular momentum J. In d space-time 
dimensions, these projectors are [9] 

= \ i 1 ^ 1 ^ - /w3rQ ) > ( 43b ) 

= \ (i^r* + 1 1 *!™) - jzi^r** ( 43c ) 

where is given in Eq. (31). For J = 0, 1, and 2, K^ v . a/3 projects the tensor A^ onto its 
trace, its antisymmetric part, and its traceless symmetric part, respectively. The factor of 
Sj(d) in Eq. (42) comes from averaging over the angular momentum states of the bb pair. 



The spin- J multiplicities in d dimensions are 

So(d) = K° map K°^ = 1, (44a) 

Si(d) = KJ^K 1 ^ = \{d - l)(d - 2), (44b) 

S*(d) = K^K 2 ^ = \{d + l)(d- 2). (44c) 



The explicit sum in Eq. (42) is over the color and spin states of the c, c, and g. In the 
expression for the amplitude A^ in Eq. (27), the only factors that depend on the spins and 
colors of the ccg are u(pi), v(p 2 ), and e°*(p 3 ). The sum over the spins and colors of the ccg 
are 

u(Pi)T a lxv( P2 K*(p 3 ) [u( Pl )T\v(p 2 )e b T *(p 3 )Y 

ccg 

N 2 - 1 

= ^ — 9arTl W 1 + m ch^fo - m ch P ] • ( 45 ) 

We have omitted terms from the sum over gluon spins that are proportional to p 3a or p 3r , 
because they give zero when they are contracted with the trace in Eq. (29) or its complex 
conjugate. After evaluating the Dirac traces in Eq. (45) and in Eq. (29), we reduce the 
contracted tensors in the differential decay rate in Eq. (42) to complicated functions of 
Lorentz scalars, which we will report later in this section. 

We wish to obtain expressions for the coefficients that are differential in the momen- 
tum of the charm quark. We must therefore integrate over the entire three-body phase space, 
except for the energy E 1 of the charm quark in the bb rest frame. The differential three- 
body phase space in the center-of-momentum frame in d = 4 — 2e space-time dimensions is 
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computed in Appendix A: 



Un) 2e r/rn „ „ dE l( lE 2 dE 3 



where A(x, y, z) = x 2 + y 2 + z 2 — 2(xy + yz + zx) and \pi\ = (E 2 — m 2 ) 1 / 2 is the magnitude 
of the three-momentum of particle i. The physical region of E ± , E 2 , and E 3 is determined 
by the delta function and by the condition — \(pl, p 2 ., pi) > 0. The physical region can be 
determined from the expression 

-\(pI,pI,P 2 3 ) = (|Pi| + |P2| + |P3|)(|Pi| + |P2|-|P3|)(|P2| + |P 3 |-|Pi|)(|P3| + |Pi|-|P2|). (47) 

We let the energies of the c, c, and g be E±, E 2 , and E 3 , respectively. It is convenient to 
introduce dimensionless energy variables Xi defined by 

Xi = Ei/E b . (48) 

We can use the delta function in Eq. (46) to integrate over x 2 . If we set E b = nib, then the 
differential phase space for ccg reduces to 

( \ 2 

^ = [(xg-r^a-co^^K]^ ^ 3 ' (49) 



where c 3 (e) is defined by 



c(e) w r (i) (2,r m 

C3(6) " r(l-e)r(|-e) " r(2-2e)' (5 ° } 



and 6*i3 is the angle between the momenta pi and p 3 : 

sm 2 e 13 = -\(pIpIpI)/{ap\pI). (51) 

The ranges of the variables X\ and x 3 are given by 

\fr < x 1 < 1, (52a) 
x 3 < x 3 < x 3 , (52b) 

where the endpoints of the x 3 integral are 

4 - - 2(1 (53, 
2-Xi^f y/xf - r 
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After integrating over the energy fractions of the c and g, we find that the differential 
annihilation rate in Eq. (42) reduces to 



8C F alA 6e 



m b 



dT[bb 1 ( i P J )^c + X] = ^^p- c^m-^tiUx^+TiXx,) dx u (54) 



where the coefficient c 3 (e) is defined in Eq. (50). The dimensionless functions f^ iv (xi) are 
defined by 

fo (r) - (d-2 + r)I (x 1 )-A[I 1 (x 1 )-I 2 (x 1 )] 

divl l) ~ (d- l)(xj -r-y ' { ' 

f i . x (d - 3)(d - 2 + r)I ( Xl ) + 4[/ 1 (x 1 ) - I 2 ( Xl )] 

LdiAXl) ~ (d-l)(d-2)(xl-rY ' 

f (d 2 - 2d - l)(d - 2 + rjlpjxj - 4(d - Wijxi) - h(xi)] 

div[Xl) ~ (d-l)(d+l)(d-2)(x 2 1 -rY ' 1 j 

where the functions I n (xi) are integrals over x 3 \ 

UXl) = I; dX \m-^^Y <56> 
These integrals, which are logarithmically infrared divergent, are evaluated analytically in 
Appendix B. They can be expressed in terms of two distributions that are singular at x\ — 1: 
the Dirac delta function 6(1 — x±) and a distribution [1/(1 — ^i)]/f that, is defined by 

/ 9{xi)[f{xi)]^dx l = I \g(x 1 )-g(l)]f(x 1 )dx 1 -g(l) f(x 1 )dx 1 (57) 

J x J x J y/r 

for any x in the interval \fr < x < 1 and any smooth function g(x\). The dimensionless 
functions fg n (xi) in Eq. (54) are defined by 

Ao / \ f 4 2(1 - x 3 )(8 + x 3 )C( Xl ,x 3 ) + 3r(4 - x 3 ) dx 3 

Ra{xi) ~ L 12(1 -x 3 y v (58a) 



n + 



ffln(zi) = A r C(xi,xs) — , (58b) 

<J Jx- X 3 

f, 2 / \ r 3+ (l-^3)(5 + ar 3 )C(a:i,X3) + 3r(2-a:3) cb 3 
fi>l) " 4 15(1 -s 3 ) 2 V (58C) 

where C(x±,x 3 ) is the function 

C(*„„)= < 1 -*> J + <? + *- 1 >' . (59) 

The results from carrying out the integrations over x 3 in Eqs. (55) and (58) are tabulated 
in Appendix C. 
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To complete the matching calculation of the coefficient dA^\ we need to evaluate the 
NRQCD matrix element on the right side of the factorization formula in Eq. (17). The bb 
states have a nonstandard normalization that corresponds to the procedure that we used 
in computing the full QCD rate gT[&&i( 3 Pj) -> c + X] [Eqs. (23) and (42)]. Application 
of that procedure in NRQCD is equivalent to the use of bb states that are normalized to 
3(2E b ) 2 /q 2 , instead of the conventional (2E b ) 2 , where q is the momentum of the b quark 
in the quarkonium rest frame. The matrix element at leading order in the nonrelativistic 
approximation is then 

(O l { z Pj)) bl ^ Pj) =m c ml (6Q) 
Substituting Eqs. (40), (54), and (60) into the factorization formula (17), we obtain 



f^O+flOrx) 



+ 



where we use 



1 2(r-l) . 
+ „L . - + log 



mi 



em 3(2 + r) 



(1 -r)A 2 



(2 + r)v^7 



(47re- 7 ) 



M € ) 
c 3 (e) 



l + 0(e 2 



5(1 -x x ) >dx 1 + 0(e), (61) 



(62) 



The explicit infrared divergence in Eq. (61) is canceled by the infrared divergence in f ^ iv (xi). 
Therefore the expression in Eq. (61) is finite at e = 0, and so we can neglect the e dependence 
in the prefactor. The only dependence on the scale A that remains appears in the bracket 
in Eq. (61). 

It is now straightforward to determine the coefficients dA^f. Our final results for the 
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differential coefficients with respect to the energy fraction %\ of the charm quark are 



dA^(A) C F a 3 s 
dx 1 



cL4j c) (A) C F al 
dxi N c 



dA^(A) C F al 
dxi 



+ 1 ( 28 



+ 



2(2 + r). ±{l-y/r)m b 1 
log — h 



9 



4 + 3r 
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log 



9 



5(1 -an) 



35xi + (2 + r) 



- + r - 3xi(l - x x ) 



1 — 

Xi 4- 



log 



Xi - a/x^ - r 



X i 2 



Xi 



Xi 



2(2 + r). 4{l-y/r)m b 1 

loe 1 

9 g v^A 18 



36 



3r , 1 
log — 



+ 



9 



5xi 



+ (2 + r) 



1 



1 — Xi 



11- 


r 


l\- 


r _ 


x\ - 


- r 



Vl-r 
5(1 - Xl ) 



1 2 - Xi + v^i - 

-^ lo §o /^^f 

<J 2 — xi — a/x| — r J 



2 

+ 9 



2(2 + r) bg 4(1 - y/¥)m b + 7 + 4r ^ ^- 



9 



40 + 21r 



73 - 89xi 



10 



+ (2 + r) 



180 
1 



1 — Xi 



log 



90 

l + y/1^7 

1 - VT^T. 



5(1 - Xl ) 



2 

xf — r 



+ -[l + r- 

1 i 2 

loe: — 

15 S 2 



Xi 



2xi (1 
' Xi + a/x| 



)] log 



Xi 



+ \fx\ 



Xi 



(63a) 



(63b) 



(63c) 



xi - \fx\- 

where [1/(1 — x\)\^ is the distribution defined in Eq. (57). 

In Eq. (63), the terms involving the [1/(1 — Xi)]y/f distribution diverge as 1/(1 — x±) as 
Xi — > 1. These singularities arises because, as X\ — > 1, the energy of the real gluon in the 
final state goes to zero, giving rise to an infrared divergence in the rate. The second term in 
the definition of the [1/(1 — xi)]^ distribution provides a negative contribution that cancels 
this divergence when one integrates over a region in xi that contains the point Xi = 1. 
Suppose that one integrates over the region x < x± < 1. Then, owing to the second term 
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in the definition of the [1/(1 — £1)]^ distribution (57), the result is dominated in the limit 
x — > 1 by a term that is proportional to — log[l/(l — x)]. Such unphysical divergences are a 
symptom of the fact that the perturbation expansion in a s breaks down in the limit X\ — > 1 
because of the appearance of large logarithms of 1 — X\. A correct treatment of the region 
near x\ — 1 would involve the resummation of logarithms of 1 — x\ (Ref. [19-24]). As x — > 1, 
real gluon emission is suppressed. Hence, the resummation of logarithms of 1 — x generally 
leads to a Sudakov factor that suppresses the rate near x — 1 (Ref. [19-24]). Consequently, 
as x — > 1, we expect the resummed distribution to turn over, rather than to diverge, and 
to approach zero smoothly at x — 1. We note that, in the rate integrated over all xi, 
logarithms of 1 — X\ do not appear, and resummation is not necessary in order to obtain a 
reliable result. 

In the limit x — > 1, the velocity expansion of NRQCD also breaks down because of 
kinematic constraints near the energy endpoint [25]. A correct treatment of this problem 
would involve the inclusion of shape functions [25] for the XbJ mesons. In general, the 
inclusion of shape functions has the effect of smearing the energy distribution near the end 
point. We expect that these smearing effects will be important for 1 — x\ less than v ~ 0.3. 
In this region, the expression in Eq. (63) should not be taken as an accurate estimate of the 
distribution. (For a discussion of these effects in the decay of the T meson into a photon 
plus light hadrons, see Ref. [26].) In the total rate integrated over x\, the velocity expansion 
is well behaved and the effects from the shape function are of higher order in v 2 . 

The resummation of logarithms of 1 — x \ and the inclusion of shape functions are beyond 
the scope of this paper. In the absence of such analyses, one should treat our results with 
caution in the region near x\ = \. 



E. Charm-quark momentum distribution 

The NRQCD factorization formula in Eq. (5) can be expressed in a form that is differential 
in the energy fraction X\ of the charm quark: 



— IXu^c + X] = dA(jC){A) I ^ { ° 8) ^\ 

dx\ dx± m\ dx\ m 2 



where the color-singlet coefficients dA j - 1 (A) /dx 1 are given in Eqs. (63) and the color-octet 
coefficient dA^ jdx\ is given in Eq. (41). 
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The momentum distribution for the charm quark can be obtained from Eq. (64) by a 
simple change of variables. It is convenient to express that momentum in terms of the 
fraction yi of the maximum momentum for a charm quark that is kinematically allowed in 
the annihilation of a bb pair at threshold: 



2 

xf — r 
1 -r ' 



The range of y± is < y± < 1. The inverse relation is 



Xi = yj(l -r)yl + r. 



The distribution in the fractional momentum y 1 can then be written as 

dT (l-r)yi dT 



(65) 



(66) 



(67) 



dyi - r )y\ + r dx l ' 

The singular distribution [1/(1 — xi)\^ in the coefficients dA ( f ) /dx 1 in Eqs. (63) can be 
transformed into a singular distribution in the variable y\ as follows. From Eq. (57) we can 
derive the identity 



1 



1 — X\ 



dx\ = < h(y 



l-2/i 



where 



%i) 



1 - y x \ dxi 



(68) 



(69) 



1-xiJ dyt 

Note that h(l) = 1. Using Eq. (66) to compute h(yi) and substituting the results into 
Eq. (68), we obtain 



1 — X\ 



dxi 



JJl 



[1 + y/(l-r)y( + 



[l + y 1 ) y /(l-r)y* + - 



l-2/i 



+ log(l + v^)5(l-yi)|dyi, 



(70) 



where the plus distribution [1/(1 — yi)]+ is defined by 

f g{yi)[f{yi)]+dyi = [\g(yi) - g{i)\f{ yi )d yi -g(i) [" f( yi )d yi (71) 

Jy Jy JO 

for any y in the interval < y < 1 and any smooth function g(yi). 

The charm-quark momentum distributions in the decays of Xbo, Xbi, and \b2 are illustrated 
in Fig. 1. For the ratio r, which is defined in Eq. (37), we choose the value r = Am 2 D /m 2 Xbj) 
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FIG. 1: Distribution of the scaled momentum y\ for the charm quark in decays of the \bj f° r J = 
(solid line), 1 (dotted line), and 2 (dashed line) for a s = 0.215, (Oi) Xb = 2.03 GeV 5 , and m b = 4.6 
GeV. 

which is equivalent in the nonrelativistic approximation that we use in the calculation, but 
more correctly reflects the physical kinematics. Here, is the average of the masses of the 
D° and D + and m XbJ is the mass of the XbJ state. 

If we use the most recent numerical values for the masses in Ref. [27], we find that this 
ratio is 0.1434, 0.1424, and 0.1419 for J = 0, 1, and 2, respectively, for the IP multiplet 
and 0.1331, 0.1325, and 0.1322 for J = 0, 1, and 2, respectively, for the 2P multiplet. 
For y < 1, the color-octet terms in Eq. (64) do not contribute at leading order in a s . 
The normalizations of the momentum distributions for y < 1 therefore depend only on the 
combination a^(Oi) Xb /mf. We choose (Oi) Xb ~ 2.03 GeV 5 , as is given in Eq. (8a). We take 
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the bottom-quark mass to be the one-loop pole mass: m b = mf° c> pa 4.6 GeV. We take a s 
to be the running coupling constant at the scale m[ polc ' ) : a s pa 0.215. The yi distributions 
for Xbo, Xbi, and Xb2 in the IP multiplet are shown in Fig. 1. As is expected from our 
discussion in Sec. HID, all three curves diverge as 1/(1 — y\) as y± — > 1. There are also 
singular distributions with support only at y± — 1 that cannot be seen in the figure. As 
we have already mentioned, the singular distributions are such that the integrals of the y\ 
distributions over an interval in yi that includes the endpoint y\ — 1 are finite. 



IV. TOTAL CHARM PRODUCTION RATE 



A. Short-distance coefficients 



The inclusive charm production rate in decays of the Xbj can be calculated by integrating 
the differential rate in Eq. (64). The integral of the color-octet coefficient in Eq. (41) is 
trivial: 



A 



(c) 



;i + r/2) v / r^ : 



(72) 



The required integrals for the color-singlet coefficients in Eq. (63) are tabulated in Ap- 
pendix D. These coefficients reduce to 
"2(2 + r) 



4° (A) 



4 C >(A) 



Cfol £ s 



4 C) (A) = ^ 



9 



■log 



1 — rjiJib 58 + 23r 



rA 



27 



5 1 + 



(73a) 



2(2 + r), 8(1 -r)m b 16 + llr 
log 



9 



rA 



27 



2(2 + r), 8(1 -r)m b 116 + 91r 
— log 



rA 



135 



r 4 l + y^7 

(73b) 

8 l + VT^T ' 

V 1 — T log . 

45 & 1 - y/T^F 

(73c) 



B. Comparison with previous results in the limit m c — >■ 

The limiting value of the color-octet coefficient in Eq. (72) as r — > is |7ra^, which 
agrees with the coefficient of n/ in the leading-order result for A s in Eq. (3c). The limiting 
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behaviors of the color-singlet coefficients as r — > are given by 

a(c), A s C F a 3 s (. 4 4, 2m fe 58 \ tf7l . 

A„ (A) — ► — — log - + - log — — , (74a) 

Ai j (A) — ► — — — log - + - log —r- • (74c) 

2 y ' N c \15 r 9 A 135 / V ; 

The coefficients and in Eqs. (74) contain logarithms of r, and they therefore 
diverge in the limit m c — > 0. In the inclusive decay rates of the x&o and the Xfe2, the 
logarithmic sensitivity of the short-distance coefficients to m c is canceled by a correction to 
the decay rate for bb — > gg from virtual cc pairs. The corrections of order o? s to the Aj from 
virtual charm quarks are given by 

^(virtual c) = _ 2 - n(0)Aj = ^ Aj lQg ^ (?5) 

where IT(A; 2 ) is the MS-subtracted quark-loop contribution to the gluon vacuum polarization 
at invariant four-momentum squared k 2 , and the coefficients A j that are nonzero at order 
a 2 are given in Eqs. (3). Then we have 

^(virtuale) = 26^^^ ^ 

N c /i 

^(virtualc) = 0> (76b) 
^(virtuale) = SC^f TUc (? } 

2 15iV c & /i v ; 

where /x is the renormalization scale associated with regularizing the ultraviolet divergence 
of the quark-loop contributions to the gluon propagator. Upon adding these terms to the 
coefficients A { ? in Eqs. (73), we see that the logarithmic dependence on m c cancels and we 
can take the limit m c — > 0. The sum of A^ and ^ virtual c ) reduces in this limit to 

Jn,(^ + ^)=^( 2ks a= t + ^-»), P7.) 

£3.(^ + ^1 =¥(^^-S' (77b) 

„3 



m. 



(c) ^ ^(virtual c) 



\ C F aj f 8 2m b 4 2m 6 116\ 

These results agree with the coefficients of n/ in the next-to- leading-order calculation of Aj 
in Ref. [9], once one takes into account the different normalization convention for (0\) Xb 
that is used in Ref. [9]. 
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C. Fraction of charm decays 



The fraction of the decays of XbJ into light hadrons that include charm is given by 



(78) 



the ratio of the NRQCD factorization formulas in Eqs. (5) and (2): 

R (c) = Af{m b ) (O^ + A^mKOs)^ 
J Aj{m b )(O x ) Xb +A,ml(0,)t b b) ' 
The short- distance coefficients in the numerator are given at leading order in a s in Eqs. (72) 
and (73). The short-distance coefficients in the denominator are given at leading order in a s 
in Eqs. (3) and (4). In Fig. 2, we show the fractions R^f as a functions of the dimensionless 
ratio ps that is defined in Eq. (11). These fractions Rj are sufficiently sensitive to p§ that 
Ps could be determined phenomenologically from measurements of the inclusive branching 
fractions of the XbJ into charm. 












0. 



2 



FIG. 2: Fractions R j of the annihilation decays of the XbJ that contain charm hadrons as functions 
of the ratio ps = mf(Os)^ b ^ / (Oi) Xb for J = (solid line), 1 (dotted line), and 2 (dashed line). 

A simple physical constraint that can be imposed on the color-octet matrix element is 
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that both the numerator and denominator in Eq. (78) should be positive. If we use the 
leading-order approximations for the coefficients, then the strongest constraint comes from 
the positivity of the numerator for J = 1. This constraint requires that p 8 > 0.032. 

V. CHARM-MESON MOMENTUM DISTRIBUTION 

In Sec. Ill, we calculated the momentum distribution of the charm quark in decays of 
the Xbo, Xbi, and xt>2- Once it is created, a charm quark will hadronize with nearly 100% 
probability into a charm hadron. The charm hadron can be the D°, D + , D s , or A c , which are 
stable under strong and electromagnetic interactions, or it can be an excited charm hadron 
whose decay products include the D°, D + , D s , or A c . The effects of hadronization make 
the momentum spectrum of the charm hadron much softer than the momentum spectrum 
of the original charm quark. 

The fragmentation of a charm quark into a charm hadron can be studied by using e + e~ 
collisions. At leading order in a s , the production of a charm hadron in e + e~ collisions 
with center-of-mass energy s/s proceeds through the creation of a c and a c with momenta 
| \J s — Ami, followed by the fragmentation of the c into the charm hadron. Since the initial 
quark has a well-defined momentum, a measurement of the momentum distribution of the 
charm hadrons provides a measurement of the fragmentation process. The CLEO and Belle 
Collaborations have measured the momentum distributions of various charm hadrons in e + e~ 
annihilation at center-of-mass energies near 10.6 GeV [28, 29]. This energy is fairly close to 
the masses of the Xb states, which are near 9.9 GeV for the IP multiplet and near 10.3 GeV 
for the 2P multiplet. The results of Refs. [28, 29] show that the effects of hadronization 
are large. It is convenient to describe them in terms of the scaled momentum y that is 
obtained by dividing the momentum by its maximum possible value. At leading order in 
a s , the distribution for the charm quark is a Dirac delta function at y — 1. The peaks of 
the distributions in y for the charm hadrons measured in Ref. [29] range from 0.59 to 0.68. 

A simple way to illustrate the effects of hadronization is to use a fragmentation approx- 
imation in which the charm-hadron momentum distribution is given by the convolution of 
the momentum distribution of the charm quark with a fragmentation function. The frag- 
mentation function D c _d{z) gives the probability distribution for a charm quark with plus 
component of momentum E\ + pi to hadronize into a charm hadron D with plus compo- 
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nent of momentum Ed + Pd = z{E\ + p±). It is convenient to scale the plus component 
of momentum by its maximum possible value in the annihilation of bb at threshold. The 
relation between the scaled plus component z 1 and the scaled three-momentum yi of the 
charm quark is 

V(l-r)^ + r + Vl-ryi 

Zl = 1 i n • ( 79 > 

l + y/l — r 

The inverse relation is 

m = ^l^fzj-r 

y 2^r^7(i + vr^)z 1 y 1 

If we neglect the difference between the mass of the quark and the mass of the charm hadron, 
there are similar relations between the scaled components zd and yp of the four-momentum 
of the charm hadron. The fragmentation approximation for the momentum distribution of 
the charm hadron can then be written as 

dT dzo f 1 dzi dyi dT 

1 — = 1 — / D{z D Zi) 

dy D dy D J z x dz l dy 1 



f j.. ^ ( Vi 1 - r )vh + r + VT^yn \ dr_ 



^{l-r)y 2 D + r J VD \ ^(1 - r)y\ + r + VI - ryi J d ^ ' 

where V(z) = zD(z). The expression for dT/dyo in Eq. (81), when integrated over y D) does 
not preserve the normalization of the total cross section J (dT/dyi) dyi, unless one takes the 
approximation of neglecting m c in comparison to m^, i.e. setting r = 0, in the relations 
(79) and (80) and in the limits of integration. In this approximation, z\ = y\ and zd = yo- 
The change in the normalization of the total cross section is negative and is of order r. This 
change is at the level of the error in the fragmentation approximation itself, which is derived 
from QCD by neglecting corrections on the order of the square of the quark mass divided 
by the hard-scattering momentum [30]. 

The Belle Collaboration determined optimal values of the parameters for analytic param- 
eterizations of the fragmentation functions for various charm hadrons by comparing their 
measured momentum distributions with the distributions predicted by Monte Carlo gener- 
ators and fragmentation functions [29]. The best fits were obtained by using fragmentation 
functions that are functions of z and the transverse momentum p±. Of the fragmentation 
functions that are functions of z only, the best fit was usually obtained by using the very 
simple Kartvelishvili fragmentation function: 

D c ^ D (z) = N D z a °(l-z). (82) 
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The fit for the D + was better than that for the D°, presumably because the momentum 
distribution for the D + has smaller contributions from the feeddown from decays of the D*° 
and the D* + . For the D + , the best fit for the exponent in Eq. (82) was a D + =4. The 



fragmentation probability for the charm hadron D. From Table X of Ref. [29], we can infer 
that the inclusive fragmentation probability for the D + , including the feeddown from decays 
of the D* + , is approximately 0.268. This fixes the normalization factor in Eq. (82) to be 
N D+ = 8.04. 

The fragmentation approximation to the D + momentum distribution that is given by 
Eq. (81) is shown in Fig. 3, where the fragmentation function is given in Eq. (82). We have 
set r = Am 2 D /m 2 Xbj , a s = a s (4.6 GeV) = 0.215, and (Oi) Xb m 2.03 GeV 5 , as is given in 
Eq. (8a), and we have chosen p 8 = 0.1. Within the fragmentation approximation, the peaks 
in the momentum distributions of the D + from decays of the XbJ are at yo = 0.53, 0.61 and 
0.58 for J = 0, 1 and 2, respectively. Also shown in Fig. 3 is the color-octet contribution 
to all three distributions, which peaks at yo = 0.79. As we have mentioned above, the 
normalization of the total cross section decreases in the fragmentation approximation by an 
amount of order r. In the present case, the fragmentation approximation decreases the total 
cross sections by about 2.6%, 1.0%, and 1.9% for J = 0, 1, and 2, respectively. 

The momentum distributions in Fig. 3 have unphysical negative values near the endpoint 
at Ud — 1- The momentum distributions very near the endpoint are dominated by the 
[1/(1 — x)]^f terms in the coefficients dA ( f > (A)/dx 1 in Eqs. (63). If we use Eq. (70) to 
transform the distribution in x\ into a distribution in y±, then the negative terms come from 
the last term in the definition of [1/(1 — y)\+ in Eq. (71). The limiting behavior as yo — > 1 
from this term in the momentum distribution is 



If we use the Kartvelishvili fragmentation function [Eq. (82)], then all other terms in dT/dyc 
vanish linearly in 1 — yp as yo — > 1- The fragmentation function D(ze>) in Eq. (83) also 
vanishes linearly in 1 — yu as yo — * 1, but the logarithm approaches — oo, and so this 
negative term dominates sufficiently close to the end point. 

As we mentioned with regard to the distributions in x\ in Sec. HID, such unphysical 
contributions arise because, near yo = 1, large logarithms of 1 — x\ cause the perturbation 



resulting fragmentation function has a peak at z — 0.8. The integral J Q dzD c ^^(z) is the 




(83) 
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FIG. 3: Distribution of the scaled momentum yu for the charm meson D + in decays of the Xbj 
for J = (solid line), 1 (dotted line), and 2 (dashed line) for p% = 0.1. Also shown is the color- 
octet contribution to the distributions (light solid line), which is the same, to within about 2 %, 
for J = 0, 1, and 2 and is already included in the other three curves. The unphysical negative 
behavior of the color-singlet contributions near the endpoint at yn = 1 might be eliminated by 
resumming logarithmic corrections to all orders, as is described in the text. 

expansion in a s to break down. We expect that resummation of these logarithms to all 
orders in perturbation theory would cure the distribution in y D of these unphysical effects. 

VI. SUMMARY 

We have used the NRQCD factorization formalism to calculate the inclusive decay rate of 
the spin-triplet bottomonium states Xbj info charm hadrons. In Eq. (5), the decay rates are 
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expressed in terms of two independent nonperturbative factors for each P-wave multiplet: 
(Oi) Xb and (Og)^ . The coefficients of these factors were calculated to leading order in 
a s using perturbative matching. Our results for the coefficients that are differential in the 
c-quark energy fraction are given in Eqs. (41) and (63). Our results for the coefficients 
integrated over the c-quark energy fraction are given in Eqs. (72) and (73). The ratios R^f 1 
of the decay rate of the Xbj into light hadrons that include charm and the decay rate into 
all light hadrons are shown in Fig. 2 as a function of the ratio p 8 of the NRQCD matrix 
elements. The ratios are sufficiently sensitive to p$ that measurements of the branching 
fraction of the XbJ into charm could be used to make a phenomenological determination of 
the (Og)x^ ■ These matrix elements could then be used to predict the partial widths into 
light hadrons for all four states in the P-wave bottomonium multiplet. 

We also calculated the momentum distribution of the charm quark from the decays of the 
Xbj- We obtained a simple approximation to the momentum distribution for charm mesons 
in Xbj decay by convolving the charm-quark momentum distribution with a fragmentation 
function for c — > D that was measured in e + e~ collisions. The charm-meson momentum 
distributions for the XbJ are shown in Fig. 3 as functions of the scaled momentum variable 
yr> for p 8 = 0.1. The CLEO-III experiment and the B factory experiments may be able to 
measure the momentum distributions of charm hadrons in Xbj decay. One unsatisfactory 
aspect of the theoretical momentum distributions in Fig. 3 is the unphysical negative be- 
havior of the distributions near the endpoint at yo — 1- We expect that this difficulty could 
be overcome by resumming logarithmic corrections to all orders in a s . The region near the 
endpoint also receives large contributions that are formally of higher order in the NRQCD 
velocity expansion. Such contributions can be resummed to all orders in v by making use of a 
shape function. The completion of these resummation calculations would allow one to make 
quantitative comparisons between theoretical predictions and experimental measurements 
of the momentum distributions of the charm hadrons that are produced in Xbj decays. 
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APPENDIX A: DIMENSION ALLY REGULARIZED THREE-BODY PHASE 
SPACE 

In d = 4 — 2e space-time dimensions, the three-body phase space is defined by 

* - - » - * - ^) ppkpPkpPk - (A1) 

where E t and Pi are the energy and four-momentum of the particle % in the final state 
with mass to*, and P = p\ + p 2 + p 3 . We evaluate d&z in the center-of-momentum frame, 
P = (y/P^,0), where the resulting expressions are most compact. In any decay with a 
three-body final state, the squared matrix element, summed over spin states, is a Lorentz 
scalar, depending only on the four momenta P, pi, p 2 , and p 3 . By using energy-momentum 
conservation, it can be seen that all possible scalar products of momenta can be expressed 
in terms of E^s. Therefore the spin-summed matrix element squared depends only on the 
energies E^. 

Integrating out p 2 and all angles except for the relative angle between p\ and p 3 , we 
obtain 

rf$3 " r(l - e )r (2 - e ) E 2 6 ^^ VP ~ El ~ E * ~ E 3 )-^-dcos0 13 , 

(A2) 

where 6*13 is the angle between pi and p 3 in the center-of-momentum frame. The angle #13 
is fixed by the energy delta function: 

E 2 = \J\Pi\ 2 + \p 3 \ 2 + 2^N cos 0i3 + m 2 2 . (A3) 
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By solving Eq. (A3) in the center-of-momentum frame, we can express sin 2 9\ 3 in terms of 
the magnitudes of three-momenta for the final-state particles: 

sm 2 9 13 = -\(plpi,pl)/(4pjpj), (A4) 

where X(x, y, z) = x 2 + y 2 + z 2 — 2{xy + yz + zx). The physical region can be determined 
from the expression 

-A(p?,P2,Ps) = (|Pi| + |Pa| + |Ps|)(|Pi| + IP2I - IPaDdPal + |Ps| - \Pi\){\Pz\ + \Pi\~ |pa|) > 0. 

(A5) 

Substituting Eq. (A4) into Eq. (A2) and changing the integration variable from cos #13 to 
E 2 , using Eq. (A3), we obtain 

In the calculations in this paper, we study the case p\ = p 2 , = m 2 c) p 3 = 0, and y/P 2 = 2E b . 

We express the energy variables in terms of dimensionless variables x« = Ei/E b . Then 

2E b (E b -E\- E 3 ) + E,E 3 (1 - Xl )/x 3 - a(x 1 ) 
cosfi 3 = i — rr\ = 77 — x > l A 'J 

IP1I-E3 K^i) 

where a(xi) and fe(xi) are defined by 

a(xi) = 1 — \x\, (A8a) 



b( Xl ) = \\Jx\-r. (A8b) 

The ranges of integrals are determined from Eq. (A5): 

y/r < X! < 1, (A9a) 

^3 < ^3 < x 3 , (A9b) 

where are defined by 

x f = l ~ X \ (A10) 

APPENDIX B: EVALUATION OF INTEGRALS I n (xi) 

In this appendix, we evaluate the integrals I n (xi) that are defined in Eq. (56): 

UX ' )= l; ^(l-orfW <B1) 
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where the bounds xf of the integral are given in Eqs. (A8) and (A10). The cosine of the 
angle #13 is expressed as a function of x\ and x 3 in Eq. (A7). By making the changes of 
variables 

I — Xi 

t = = a(x{) + b(x 1 )cos9 13 (B2) 

^3 

and using the relation for cos 6*13 in Eq. (A7), we can parametrize the I n (xi) as 



{i-xi) 1+2£ Ja( X1 )- b ( Xl) (1- cos 2 e 13 y 

b(xi) f 1 [a(xi) + b(xi)x] n+2e dx 



(B3) 



In the second line of Eq. (B3), we used Eq. (A7). It is evident that the t or x integrals in 
Eq. (B3) are finite. The divergent part is contained in the factor 1/(1 — x 1 ) 1+2e , which is 
manifestly logarithmically divergent in the infrared limit X\ — > 1. The integral of that factor 
over xi is proportional to — l/(2e). The evaluation of the integrals I n (xi), keeping the full 
e dependence, is quite involved. However, in order to compute the pole in e and the finite 
term, we need only to expand the coefficient of 1/(1 — x\) 1+2e in Eq. (B3) to order e: 



{i-x l y+^ {j a{xi) - b{xi) 



-eb(xi) j [a(x x ) + b(x l )x\ n \og(l - x 2 )dx^ + 0(e) 

[a(xi) + 6(xi)] n+1+2e - [a(x 1 ) - fe(xi)]" +1+2e 



(l- Xi y+^{ n + l + 2e 2 " n(Xl) 

+0(e), (B4) 

where the i n (xi) for n = 0, 1, and 2 are given by 

i ( Xl ) = -26(a;i)(l-log2), (B5a) 
ii(xi) = -2a(xi)6(a:i)(l-log2), (B5b) 
i 2 ( Xl ) = -26(x 1 ){i[^i)] 2 (4-31og2) + [a(x 1 )] 2 (l-log2)}. (B5c) 



It is convenient to rewrite the divergent integral as a linear combination of finite integrals 
and a singular integral involving a delta function. As we have noted, all the factors except 
for 1/(1 — xi) 1+2e are regular functions of x\. We denote the factor that is the coefficient of 
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1/(1 — xi) 1+2t by f(xi). Therefore, we wish to study the integral 

f{x 1 )dx 1 



\l+2e ' 



(B6) 



(1 -x Xt 

where f{x\) is regular for any X\ G [y/r, 1]. One can separate the divergent contributions to 
the integral from the finite piece as follows. 



1 



**! + /(!) f TT 



|l+2e' 



(B7) 



, (l-xi)^ 

The first term on the right side of Eq. (B7) is finite in the limit e — > 0. The second term on 
the right side of Eq. (B7) is singular in the limit e — > 0: 



r — 



^(l-an)^ 2e(l- V^) 2e ' 
Substituting Eq. (B8) into Eq. (B7), we obtain 

-f( Xl )-f(l) f( Xl )5(l - Xl ) 



(B8) 



I = 



dxi 



V7 



(1 - Xl y+ 2 * 2e(l - V^) 2e 



(B9) 



Expanding the e dependence in the first term in Eq. (B9), we have 



5(1 -x ± ) 



(l-x 1 ) 1 + 2 " 2e(l-y/r) 



2c 



+ 



1 — Xi 



2e 



log(l - Xi) 
1 — Xi 



+ 0(e 2 ), (BIO) 



where the distribution [1/(1 — ^i)]^ is defined by Eq. (57). Retaining the first two terms 
in Eq. (BIO), we obtain 



(xj — r-y 

h(xi) 
(x\ — r) e 

h(xi) 



-l + log2 ) sA ~ r I U(, ) 



5(l-xi) + 



1 — Xi 



\Jx\-r + 0(e), 

(Blla) 



5(1 -xi) 



+ 



1 


1 




2 


1 — Xi 


J 



xj -r + 0(e), 



(Bllb) 



1 1 



+ 



24e 12 
2 — xi 4 — r 



1 r 



— + — log2 (4-r)-- + — 



4 12 



v / T^ + L 2 (r)jcKl-x 1 ) 



1 — 



(Bile) 



34 



The functions L n (r) are given by 

L n (r) = \- [(a + b) n+1 log(a + b) - (a - b) n+1 log(a - b)] 

+^Vl) [(a + &) n+1 - (« - log [(1 - V~rf{l - r)] , (B12) 

where a ± 6 = a(l) ± 6(1) = |(1 ± vl — r). These functions vanish in the massless limit 
r — > 0: L n (0) = 0. The explicit expressions for Lo(r), -^i(r), and ^(r) are 



EIZL 4(1- y^) 2 (l-r) 1 l + y^7 
L (r) = log — - — — log , B13a 



rf\ y_[EL\ 4(1- y^) 2 (l-r) 2-r 1 + y/T^T 

L i( r ) = 5 lo S 5~ lo Si 7^=' (B13b) 



4 & r 8 & l-v^ 



r 



(4-rVTHZ, 4(1- y^) 2 (l-r) 4-3r 1 + VT^ 

L 2 (r) = ; log — — 1 log ; B13c 

y ' 24 6 r 24 6 1 - v ; 

APPENDIX C: EVALUATION OF INTEGRALS OVER x 3 

In this appendix, we report the results of carrying out the integrations over X3 in the 
components of the bb differential widths in Eqs. (55) and (58). 

The integrals I n (xi) that appear in the infrared-divergent functions f ^ iv (xi) in Eqs. (55) 
are defined by integrals over x 3 that are evaluated in Appendix B. Making use of these 
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results, we find that 

fdiv(zi) = 



2(2 + r) 
9 



-- + log2 



5 + r\ ; — 



2 + r 



L (r) - ^(r) - L 2 (r)}^ 5(1 - Xl ) 



+- { l-2xi + (2 + r) 



1 




1 — Xi 





2(2 + r) 
9 



+ 



2e 
2 + r 



+ log2 



xf — r, 



7-4r\ 

+ —rr- I vi - 



(Cla) 
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6 Lo(r) + -[L 1 (r)-L 2 (r)]}S(l-x 1 ) 



2 
+ 9 



1 - 2xi 



+ (2 + r) 



1 — Xi 



x\ — r, 



(Clb) 



2(2 + r) 
9 



41 - 8r \ ^ — 
270 1 



+^5^A)(r) - ^[Li(r) - La(r)] |> 5(1 - x : ) 



+ 



2 J 1 - 2xi 
10 



+ (2 + r) 



1 




1 — Xi 





xf — r, 



(Clc) 



where the distribution [1/(1 — Xi)]^ is defined in Eq. (71). 

The infrared-finite functions rg n (xi) in Eqs. (58) are defined by integrals over X3 that are 
straightforward to evaluate. The results of carrying out these integrations are 

2y/x\ - r 



r§n(*i) = 



3 

-7 log 



(9- llxi) + 



- + r - 3xi (1 - xi) 



log 



xi + y/xj-r 



Xi 



y/xj ~ r 



1 2 - xi + y/xj - r 



rL(^i) 
fL(^i) 



6 

X 1 y/xf 



2 — xi — a/x^ — r 
2 - xi 



--lie". 

3 3 2 — xx — \fx\ — r 



(C2a) 
(C2b) 



15 



■(24 - 29xi) + - [1 + r - 2xi(l - Xi)] log 



Xi 



Xi 



1 2 - xi + y/xj-r 



15 



2 - Xi - y/xf 



(C2c) 



36 



APPENDIX D: EVALUATION OF INTEGRALS OVER x x 



In this appendix, we tabulate the integrals that are required to obtain the inclusive short- 
distance coefficients A^f in Eq. (73) from the short-distance coefficients dA^f in Eq. (63) 
that are differential in x\. 

Some of the basic integrals over x\ are 

d.r i ./• \ / .ri — r = . 

VT^7 



I dx x x x sjx\-r = Ul-rf' 2 , 



/ dx\ \l x\ — r 



r 1 + \J\ — r 



(Dla) 
(Dlb) 



/ dx\ J x\ — r 



1 — X\ 



I + ilog^-log(I + v^) 



1 , 1 + 
— log , 

2 s i_ v /r^7 



(Die) 



The integrals over x\ of the infrared-divergent functions f^ v (xi) given in Eq. (CI) are 



dx^^Xi) 



^lfL^i) 



2(2 + r) 




1 


9 




"27 


27 


— r 


2(2 + r) 




1 


9 






2(4 + 5r) 


Vi 


27 




2(2 + r) 




1 


9 




"27 


8(5 + 4 


r) 


Vi 


135 





log- 



8 + 3r 
' 18 

r 

- + 3>c 



log 



'l-r 






r 




r 


1 — r 


Vl -r 


1 + y/l-r 


1 - VI - r' 







(D2a) 



(D2b) 



VT 



40 + 21r , 1 + Vl^r 
log 



90 ~° 1-^/1=7' 

The integrals over xi of the infrared-finite functions f g n (xi) given in Eq. (C2) are 

! 8(2 + r) /- 3 + r 1 + VT^7 

VI - r + 



(D2c) 



/ ^lffi n (Xi) 



9 

2 + r 



3 log l-VT^7' 
1 + VT^7 



-i 



9 Vr37+ 6 K i-VT37' 



/ rfxif| n (xi) 



22 + 23r , 8 + 7r, 1 + VT^ 

■VI — r H ^— log ■ 



45 



30 



1 - VT^ 7 



(D3a) 
(D3b) 
(D3c) 
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